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We derive a general and simple expression for the time-dependence of the position operator of 
, a multi-band Hamiltonian with arbitrary matrix elements depending only on the momentum of 

(N . the quasi-particle. Our result shows that in such systems the Zitterbewegung like term related to 

f— I ' a trembling motion of the quasi-particle, always appears in the position operator. Moreover, the 

, Zitterbewegung is, in general, a multi-frequency oscillatory motion of the quasi-particle. We derive 

a few different expressions for the amplitude of the oscillatory motion including that related to the 
' Berry connection matrix. We present several examples to demonstrate how general and versatile 

, our result is. 
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, Schrodinger in his original paper has predicted a 'trembling' or in other words a rapid oscillatory motion of the 
^ • center of the free wave packet for relativistic electron However, the Zitterbewegung is not strictly a relativistic 
D ; effect but can be observed in spintronic systems as well [7|. This work has initiated many other works |8l-l32l| 

with an aim to demonstrate the appearance of the Zitterbewegung not only for relativistic Dirac electrons. In these 
works a common feature is that the oscillatory motion of the free particles can be described only by one frequency. 
Cd ' In our previous work we showed that for a wide class of Hamiltonians related to for example spintronic systems and 
graphene, the Zitterbewegung can be treated in a unified way i33i] . Here the basic idea was that the Hamiltonian of 
several systems can be mapped to that modelling the precession of a virtual spin in an effective magnetic field. The 
coupled equations for this virtual spin precession and the orbital motion of the quasi-particle can easily be solved. 
Thus, the Zitterbewegung is arising because the virtual spin and the orbital motion for the quasi-particle are coupled. 
^ , Our work suggests as natural question whether the phenomenon of the Zitterbewegung also arises for an even more 
• general Hamiltonian. 

' In the present work we extend the Zitterbewegung phenomena to a broader class of quantum Hamiltonians for free 
^ (quasi-) particles. In particularly, we derive a general and simple expression for the time-dependence of the position 
operator x(i) for a multi-band Hamiltonian given by 



o 



H 



i?2i(p) if22(p) ... H2n{p) 



Vi?nl(p) i/„2(p) ... Hnnip)J 



(1) 



where each matrix element is a differentiable function of the momentum p of the particle itself and n > 2 is the 
number of degrees of freedom of the system. From our general expression for the position operator x(f) we shall show 
that i) the Zitterbewegung always appears for systems given by the Hamiltonian ([T]) , i) for n > 2 the Zitterbewegung 
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X . . . 

^ , is in fact a multi- component oscillatory motion of the free quasi-particle, ii) for n = 2 we recover the results obtained 
■ - - I earlier in the above mentioned references. 

To find the time dependence of the position operator x(i) of the quasi-particle in Heisenberg picture one needs to 
calculate 

x(i) = e*^*x(0)e-^^*, (2) 

where x(0) is the position operator at t = 0, ie, it equals to the position operator in Schrodinger picture. Because 
the momentum operator p is a constant of motion we now work in the subspace of the Hilbert space for which 
the momentum p is fixed. Calculating the right hand side of equation ([2]) the crucial step is to decompose the 
Hamiltonian ([T]) into a sum of projection operators: H — J^a ^aQa, where Ea is the ath eigenvalue of the Hamilton 
operator at a given momentum p, and Qa are projection operators satisfying the following relations: QaQt — ^ab Qa 
and X)a ~ where /„ is the nx n unit matrix. The position operator at time i = in Schrodinger picture 
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and in momentum representation is x(0) = i^^- Consider the operator U = e r ^* which is only a function of the 
momentum operator p. Then equation ^ can be rewritten as 



x(t) U-^x{Q)U = U-^ [x(0),t/] + t/-iC/x(0) 
dU 
dp 



871 

= x(0)+iftC/-i^, (3) 



where we have made use the relation [x(0),F(p)] = ih^^^. Decomposition of the Hamiltonian ([ij into a sum 

of projection operators makes possible to write that e^^^* — J2a^^^ ^"^Qa- Now, substituting these operators 
into equation ([3]) and using the orthogonality relations QaQt = ^ab Qa it yields the time dependence of the position 
operator: 

X(0 = ^{0)+Y,'^aa+tY,^aQa + 

a a 

+ ^h*^^^ (4a) 

a b^a 

Va(p) = , Zab[p) ^thQa-^, (4b) 

and oJab = are the so-called beating frequencies. Here we call Vq as partial velocities and Z^fc as Zitterbewegung 

amplitudes. This is our central result in this work. 

The interpretation of the different terms in is as follows. The first term is the initial position of the quasi-particle. 
In contrast to the usual dynamics (for systems with one degree of freedom), the second and the fourth term are entirely 
new. The second term is a displacement of the position operator independent of time. The third term describes the 
motion of the quasi-particle with constant velocity which is, in general, not equal to any of the partial velocities 
Va. Finally, the Zitterbewegung stems from the last, oscillatory term which describes the oscillatory motion of the 
quasi-particle. The phenomenon of Zitterbewegung is similar to the beating effect with different frequencies in the 
classical wave mechanics. The Zitterbewegung is a direct consequence of the coupling of different energy eigenstates 
for systems with more than one degree of freedom. These terms in x(i) are inherent of the Zitterbewegung and are 
expressed via the projection operator Qa related to the given Hamiltonian. 

Equation (0]) is the most general form for describing the phenomenon of the Zitterbewegung. Our result shows 
explicitly that the oscillatory motion (the last term in equation (j4ap ) is a superposition of individual oscillatory 
motions with frequencies corresponding to all possible differences of the energy eigenvalues of the Hamiltonian ([1} . 
Thus in the most general case the Zitterbewegung describes a multi-frequency oscillatory motion of the quasi-particle. 
This multi-frequency behavior of the Zitterbewegung has first been shown by Winkler et al. in Ref. for two specific 
systems, namely for the Kane modell and Landau- Rashba Hamiltonian. However, the most clear manifestation of this 
multi-frequency behavior of the Zitterbewegung for the general Hamiltonian ([T]) can be seen only in our main result 
(HI). In general, the Zitterbewegung cannot be described by only one frequency (this is the case only for systems with 
two different eigenenergies) but all of the differences between the different energy eigenvalues corresponding to the 
beating frequencies appear in the time dependence of the position operator. 

Sometimes in the explicit calculation of the position operator x(i) it is more useful to use a different form for the 
Zitterbewegung amplitudes Zab given by Taking the derivative of the Hamilton operator H ~ E^Qc and the 
orthogonality relation QcQb — Scb Qb with respect to the momentum p one can easily show that (for details see Sec. A 
in the Appendix) 

Zab^tfl ^ , (5) 

Hjb — t^a 

valid for a ^ b. Thus in the calculation of the Zitterbewegung amplitudes instead of knowing the derivative of the 
projection operators with respect to the momentum one needs to take only the derivative of the Hamiltonian. Another 
form of the position operator x(t) is given in the Appendix. 

We now consider several examples to demonstrate how versatile our result is to study different systems known in 
the literature (for more details see the Appendix) . Regarding the Zitterbewegung most of the systems studied in 
the literature are described by a Hamiltonian with only two different eigenvalues |lM33|. In such systems either the 
Hamiltonian itself is a 2 x 2 matrix or the dimension of the Hilbert space is more than 2 but the eigenvalues are 
degenerate and the Hamiltonian has only two different eigenvalues. Thus for such systems it is useful to derive an 
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alternative form for the time dependence of the position operator given by Eq. (|4]). Now the Hamiltonian in terms of 
projectors reads H = E+Q+ + where Q± are the projection operators satisfying the usual relations mentioned 

above and E± are the two eigenvalues of H. Introducing the operator T = — Q_ it is obvious that = I. Note 
that in the mathematical literature the operators satisfying this relation are called involutary operator related to the 
mirror image in geometry. The Hamiltonian can be rewritten as H = e I -\- {huj/2) T, where e = [E^ + _E_)/2 and 
Lu = {E^ — E-)/h. Moreover, it is clear that Q± = (/± T)/2. Then using Eq. Q one can easily show that 



x(i) 
W 



= x(0)+Wt + Z(t), where 



de ^ 1 dhoj ^ 

— / H T 

dp 2 dp ' 



lit) = 2^mM- + -(l 



dT 

cos{ujt))T—. 

op 



(6a) 
(6b) 

(6c) 



From this result it is clear that in the oscillatory part of x(t) there is only one frequency component. 

Equation ([6]) can easily be applied to the original Schrodinger's Zitterbewegung and we find the same results as 
that by Schrodinger (see Sec. B of the Appendix). Another example is the Luttinger Hamiltonian 3J, |35| given by 



H 



1 

2™ 



71 + 2 T2 



p2- 272(pS)^ 



(7) 



where p = (pj:,Py,Pz) is the vector of the momentum operators, S — (Sx, Sy, Sz) represents the spin operator with 
spin 3/2, m and 71.2 are parameters of the model. Using Eq. ([5]) the position operator for Luttinger Hamiltonian can 
easily be derived (for more details see Sec. C of the Appendix): 



^n^^ hl±^2lT 272 (pS)M , ^ ■ . ,J P(v^ S(pS) + (pS)S ^ 

x t = X + ^ — I4 — pt + sm{ujt) J 

\ m m / V P 2 p^ / 

. (P X S) (pS)^ + 2 (pS) (p X S) (pS) + (pS)^ (p X S) 

where cj = E^ — E^ = [2^2 I'm) p^ . Note that this result agrees with that obtained by Winkler et al. in Ref. [l^ and 
by J. Schliemann in a private communication using a direct calculation of the right hand side of equation 

We also consider a non-trivial example for the Zitterbewegung not known in the literature, namely the Zitterbe- 
wegung for bilayer graphene. The Hamiltonian for bilayer graphene in the four by four representation is given in 



Refs. |36l |37|. Including the trigonal warping j37H39l | the position operator x(t) is more cumbersome but its structure 
and the steps of the derivation are similar to the case when the trigonal warping is omitted. Therefore, we now ne- 
glect the trigonal warping. The position operator x(t) can be derived using Eq. ^ but to obtain the Zitterbewegung 
amplitudes it is more effective to use Eq. ([S]). The results is quite lengthy thus here we only refer to Sec. E of the 
Appendix. This is a non-trivial example for the Zitterbewegung. Since the Hamilton operator for bilayer graphene 
has four different eigenvalues, we have six values of the energy differences. However out of these six values there are 
only four different ones. Therefore, the number of beating frequencies is only four. In this example it is clear that the 
oscillatory motion of the electron is a superposition of individual oscillatory motions with four different frequencies. 

Recently, for specific systems the connection between the Zitterbewegung and the Berry phase has been noticed 
and investigated by Vaishnav and Clark [23], and Englman and Vertesi (28!|. We now show that the oscillatory terms, 
ie, the Zitterbewegung amplitudes in the position operator have a close relation to the Berry connection matrix 
appearing in the expression of the well-known Berry phase (40| even for a general Hamiltonian ([l} . To this end we 
present another form for position operator in terms of the eigenvectors of the Hamiltonian. 

The projection operator can be expressed via the eigenvectors \ua.s{p)) of the Hamiltonian operator: Qa{p) ~ 
X^s \'^a.s{p)) {ua.s{p)\, whcre s denotes the different eigenvectors in a subspace with the same energy eigenvalue Ea- 
Then equation ^ can be rewritten as 

x(t) = x(0)+i^ Vfc|Mfc(p))(iife(p)| 



k 



-1) Afc,(p)|ufc(p))(uKp)|, (9a) 



k,l 



Afci(p) = ih{uk{p)\-^\ui{p)). (9b) 
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Here A^; is the so-called Berry connection matrix. The index k labels the eigenvectors of the Hamiltonian with taking 
into account their multiplicity. 

For systems with precessing spin in an effective magnetic field it turns out that to study the Zitterbewegung Eq. Q 
is more appropriate than Eq. This is demonstrated in Sec. D of the Appendix, where we find the same result for 
the position operator as that we derived before using a different approach [33j . 

So far we concentrate on the structure of the Zitterbewegung for general Hamiltonian ([T]). However, the observation 
of the Zitterbewegung experimentally is more difficult problem. It is well-known that the spatial size of the trembling 
motion of the relativistic electron predicted by Schrodinger is of the order of Compton wavelength, and its frequency 
is far beyond the present experimental possibilities [H. The experimental observation of the Zitterbewegung in the 
non-relativistic quantum regime such as in semiconductors with spin-orbit couplings is much more promising. For 
example, Vaishnav and Clark psj . and Merkl et al. [3l| have proposed an experiment for observing Zitterbewegung 
using ultra cold atoms, while Rusin and Zawadzki have proposed an experiment for observing Zitterbewegung probed 
by femtosecond laser pulses in graphene [3^ . Very recently, Gerritsma et al. have performed a quantum simulation 
of the Dirac equation using a single trapped ion and observed the Zitterbewegung 41 1 . 

One of the difficulty of observing the Zitterbewegung is the lack of time resolved probes. The initial state, in general, 
is a superposition of different momentum-eigenstates: \'^o) = J cfip X]fc'^'=(P) Wk{p))- Therefore, the expectation 
value is x{t) — (^Po | x(t) | 5'o) = / d^p J2k i "^fcCP) ^((p) (Mfe(p) I |'«i(p)) which involves the integration over the 
momentum p. Since in Eq. (j4]) the beating frequencies uJki{p) depend on the momentum p the integration over 
the momentum p in x(i) may result in a strong suppression of the Zitterbewegung in time. This problem can be 
circumvented if at least one beating frequency is independent of the momentum p of the quasi-particle. This is the 
case, for example, for bilayer graphene as shown in Sec. E of the Appendix. One beating frequency is constant and in 
the detectable regime uj = ji/h ^ 0.6 fs~^, where 71 ~ 0.4 eV is the strongest interlayer coupling between two carbon 
atoms that are on the top of each other 36, 37, 37-3^. The amplitudes of the trembling motion will be investigated 
in the near future. Our main aim in this paper is to establish a general theory for the Zitterbewegung. On the 
other hand our general theory can be a good starting point to search for systems that are realistic for experimental 
observation of the Zitterbewegung. 

Conclusions. — We presented a general theory for Zitterbewegung and derived a general and simple expression for 
the position operator x(i) in Heisenberg picture and in momentum representation, and for a given system it can easily 
be calculated. In contrast to systems studied in the literature Il9l - l33l | the Zitterbewegung is a universal 

phenomenon and it always appears in the quantum dynamics of a system of quasi-particle with more than one degree 
of freedom. Our main result ^ shows that the Zitterbewegung, in general, is a multi-frequency beating effect in 
Heisenberg picture and has a close relation to the Berry connection. We believe that our work presented here provides 
a better understanding and experimental guide for the Zitterbewegung studied intensively in the literature. 
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APPENDIX 



We now apply our main result given in the main text for Zitterbewegung to several systems known in the literature 
and also calculate the operator x(i) for non-trivial systems previously not studied in the context of Zitterbewegung. 
Below we also present the details of the derivation of the results mentioned in the main text. Our examples show how 
powerful our general expressions (H]) or ([S]) are for studying the phenomenon of Zitterwbewegung. 



A. Other equivalent forms for the position operator 

In this section we present the derivation of the result given by Eq. (O for the position operator x(t). Using the 
relation ih Qa — following from J^b Qb = In, Eq. (|3]) can easily be rewritten as 

x(i) - x(0) = tJ2^Q, + ^hJ2 (e^""^ * " l) Qa (10) 

We now present another form for the Zitterbewegung amplitudes Zab given by Eq. (U) . Taking the derivative of the 
Hamilton operator H = EcQc with respect to the momentum p we have 

Similarly, from the orthogonality relation QcQb — Scb Qb we find that 

dQc ^ , n ^*56 . dQb , . 

~Q^^b + Qc-Q^ = Ocb-^- (12) 

Thus using Eqs. (jlip and ([T2|) for a given a and b we may write: 

n n \-f^^^nnnj.T?n^Q^n\ n n ( n 

Wa — > I -T{ — HaHcHb + ^a—^ ^b — — VaVb + > , i^cHa Ocb ^c—^ 

ap ^ V c'P ap / op ^ \ op ap 

= <5ab^Qa + (£;b-K)ga^, (13) 
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where we have made used of the orthogonahty relations QaQb = ^ab Qa- Hence we find that for a b 

dQb _ Qai§Qb 

Now from this relation it is simple to obtain the alternative expression for the Zitterbewegung amplitudes given by 
Eq. ©. 

To get deeper insight into the phenomenon of the Zitterbewegung it is useful to present another form for the time 
dependent operator x(t). For the given index a and b in Eq. (jlOp we now introduce the operator Tab = Qa — Qb- Then 
it follows that T^^^ = Qa + Qb and T^^^ = Tab- The operators satisfying the latter equation are called weak involutary 
operators. It is obvious that an involutary operator is always a weak involutary operator. Note also that Ti,a — --Tab- 
From the definition of the operator Tab we find that Qa — \ (T^^, + Tab) and Qb = ^ ^ab ~ Tab)- We now express the 
second and third terms in Eq. (|10p in terms of the operator Tat- Then one can easily show that 



dhuja. 



op op an 

a.b a 



(15) 



^ ^*5b , ^ dQa 1 2 9Tab , . 

^ dQb ^ dQa 1^ dTab^ 

— n^ab—r^ J^ab, [i- < ) 

op Op 2 op 

where /„ is the n x n unit matrix. Then it is easy to rewrite Eq. (jlOp as 

ij% dT 

x(i) = x(0) + Wt+ — [i sin(lUa6 t) + Tab COs{uJab t)]Tab^Tab, (18a) 

a b^a 



= x(0) +Wt + -y (e*"-^-* _ 1) ^ Tab, where (18b) 
4 ^ ' ap 

a, 6 

n Op zn ^-^ Op 

a.b 

In our general case, when there are more than one frequency component in the oscillatory term owing to the fact 
that the system has more than two different energy levels, the time dependent form of the oscillatory term has the 
same exponential form between the corresponding energy levels as in Schrodinger's case (see the next section of the 
Appendix). In summary, Eqs. (fTO|) and ([T8| are the two alternative forms of our main result ^ for the position 
operator. 



B. Revisiting Schrodinger's original derivation 

Our expression given by Eq. ([6]) can directly be used to study the Zitterbewegung for free electron described by the 
Dirac equation. This problem was originally examined by Schrodinger [ij. The Dirac Hamiltonian reads 

H = cap + mc^l3, (19) 

where a and /? are the usual 4x4 matrices satisfying the relations of — — ctj = 0^ = 1, ctiP + Pea — and 

a^afc + akOii — 2Sik- One can show that = i?^, where E — \J c^p^ + (mc^)^. Introducing the operator T = H/E 
it is clear that T is an involutory operator since T^ = TP' jE^ = /, where / is a four by four unit matrix. Then the 
Hamiltonian can be written as H — ^T. Thus, from Eq. ^ the time dependence of the position operator is 

ft duj „ h , ^ dT h , , ^ , „dT 
x(<) = x(0) + -—Tt+- smiu;t) ^ + ^ (cos(^t) ^ 1) . (20) 

We need to evaluate the gradient of uj and T with respect to p: 

dT d{H/E) v-V 

— = ^-^ = — — , where v = ca. 22 
op op E 



7 



Here V is the classical relativistic velocity of the particle with momentum p and energy H. Then x(<) becomes 

x(t) = x(0) + V t + ^ sin(^i) + ^ (cos(^i) - 1) T 

= x(0)+Vi + *ft(v-V)^^^^. (23) 

This results agrees with that given by Schrodinger in his original work on the Zitterbewegung of the relativistic 
free electron. 



C. Luttinger Hamiltonian 



We now present a non-trivial example for calculating the Zitterbewegung in case of Luttinger Hamiltonian given 



by 



H = 



2m 



7i + §72)p'- 272(pS)^ 



(24) 



where p — {px,Py,Pz) is the vector of the momentum operators and S = {Sx,Sy,Sz) represents the spin operator 
with spin 3/2, while m and 71^2 are parameters of the model 34, 3^ (in this section we take h = 1). The Hamiltonian 
can be expressed in terms of the projection operators Ql and Qh as 



Ql(p) 

Qh{p) 



H - Eh{v)Qh{\>) + El{p)Ql{i>), where 

^-J - — 
8 ^ 2p2 



(25) 
(26) 
(27) 



and I4 is the 4x4 unit matrix, and the double degenerate eigenvalues are £'l(p) — ''^tm^ ^^^'^ ^nip) = ''^2m^ 
corresponding to the light-hole (L) and the heavy-hole (H) bands. The projection operators Ql and Qh satisfy the 
usual orthogonality relations. 

To find the time dependence of the position operator we use Eq. ([6]) valid for Hamiltonian with two eigenvalues. Take 



Ql = Q+ and Qh - Q- then we find T = 0+ - Q_ = | /4 - ^ 
Furthermore, we need the derivative of T with respect to p: 



E++E^ 



dT _ 2p(pS)-' _ S (pS) + (pS) S 
dp p^ p2 

Now the time dependence of the position operator reads 

x(t) =x(0)+W< + Z(t), 
where W and Z{t) can easily be obtained from Eq. ([6]): 



^p2anda. = i?+-£;_ = ^p2. 

(28) 

(29) 



W 



71 



:72 



1 



p/4 

dT 



272 p(psr 

m p2 

i 



Z{t) = ^ sin (ut) (cos(wi) - 1) (pS)^S (pS) - (pS) S(pS)^ + (pS)^S - S{pSf 

3sion Z(t) can be further simplified 
the commutation relations [Sj, Sk] ~ iSjuSi) 



(30) 
(31) 



2 ap 4p4 

The expression Z(t) can be further simplified using the following identity for the spin operators (can be proven using 



S(pS)-(pS)S = ipx S. 
Finally, the position operator (j29p for Luttinger Hamiltonian can be written as 



(32) 



x(t) = x(0) + 



71 



§72 ^ 272 (pS)^ 
J-i 



p i + sin {iot) 



p(pS)' S(pS)-K(pS)S 



2p2 



-|-(1 - cos(cjt)) 



(p X S) (pS)" + 2 (pS) (p X S) (pS) + (pS)^ (p X S) 
4p4 



(33) 



Note that this result agrees with that obtained by Winkler et al. in Ref. 19], and by J. Schliemann in a private 
communication using a direct calculation of the right hand side of Eq. ^ . 
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D. Spin in an effective magnetic field 



In this section we consider the Zitterbewegung for systems mapped to a system of virtual spin in an effective 
magnetic field. Such classes of Hamiltonian have been previously studied by the present authors in Ref. [s^ using a 
different approach. For many systems 

0, 0-13 [iS-Iii 

the Hamiltonian can be written in a quite general form: 



H = e(p)/ + f2(p)S, (34) 

where / is the unit matrix in spin space, and the system is characterized by the one-particle energy dispersion e(p) 
and the effective magnetic field r2(p) coupled to a virtual spin S with ma gnit ude S. Here we assume that e(p) and 
r2(p) are differentiable functions of the momentum p = {px,Py,Pz)- In Ref. [33| we listed a few systems (together with 
the effective magnetic field J^(p)) that are currently intensely studied in spintronics, and in the research of graphene 
and superconductors. To find operator x{t) for 5 = 1/2 one can use our result given by Eq. however for S ^ 1/2 
the derivation is more subtle. We now highlight the main steps in this calculation for the case oi S 1/2. 

Now it is more appropriate to start from the alternative form for x(t) given in the main text by Eq. ^ that involves 
the Berry connection matrix. Solving the eigenvalue problem for Hamiltonian (j34p we find that the eigen-energies 
are Em = e{p) + h^m, where fl = Vi}F and m = —S, —S + 1, . . . , S — 1,S labels the eigenvectors |m) of H which 
are the same as the eigenvectors of operator Sz , the component of the spin operator S with spin quantization axis z 
pointing along vector ft. In what follows, it is useful to introduce the unit vector n = 0/J7. Then (nSo)|m) — m\m), 
where So is the spin operator S in Schrodinger picture, ie, it is time independent. The Hamilton operator becomes 
H = Em=-s ^™ \m){m\. 

In the second term of Eq. (I^a)) the velocity operator W \uk){uk\ can be obtained as 

^ = 1^ ^'""^^""'"9^ ^ |TO)(m| + ^ 2^ TO|m)(m| = — /+— -(nSo), (35) 

m— — S m—~S ra— — S 

where we have made use of the fact that the eigenvectors |m) form a complete set. 

Similarly, we can calculate the Berry connection matrix Am„ in Eq. (j9bp . Using the general relation 

d f- \n) 
(m| — |n) = — 2-^, m^m', (36) 

op J^n — tLm 



(^.L„ = i? — E K,i{m\Si\n), (37) 



the mn matrix element of the jth component of the vector operator A can be expressed as 

ih \ ^ 

valid for m ^ n and where the matrix K is defined as Kjk = dflk/dpj. Then employing the well-known relations 
between the spin operators , Sy and , one can evaluate the matrix elements of the spin operator Sj in the above 
equation. After some algebra the time dependence of the position operator becomes 

x(i) = x(0)+Wt + Z(t). where (38) 

Z(t) = 5il^K(/-non)So+ '~™"' K(.ixSo). (40) 

and n o n denotes the outer or direct product, ie, (n o n)^-, = njUk- Note that this result agrees exactly with that 
derived with a different method in our previous publication [33j. 



E. Bilayer graphene without trigonal warping 

We now present a non-trivial example not known in the literature, namely the bilayer graphene system. In this 
case we show that the oscillatory motion of the electron is a superposition of individual oscillatory motions with four 
different frequencies. 
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The Hamiltonian for bilayer graphene in the four by four representation reads 36l. |37| 



H 



















vp+ 













71 


\vp+ 





71 


/ 



(41) 



where p± = p^zt ipy (note that p+p-. — p'^+Py — p^). To have a result independent from the choice of the coordinate 
systems it is convenient to introduce the vector p — {px,Py, 0). Here v — \/3ajo/{2h), where a is the lattice constant in 
the honeycomb lattice, 70 is the intralayer coupling between nearest neighbor carbon atoms, while 71 is the strongest 
interlayer coupling between two carbon atoms that are on the top of each other. In the above Hamiltonian we neglect 
the trigonal warping [s^-H^ since the other interlayer couplings are much smaller than 71. We have also calculated 
the time dependence of the position operator when the trigonal warping term is included. In this case the result is 
more cumbersome (not presented here) but the structure of the position operator and the steps of its derivation are 
similar to that presented below when the trigonal warping is omitted. Therefore, for the sake of simplicity we now 
only consider the bilayer graphene without trigonal warping. 

In what follows it is convenient to work with dimensionless Hamiltonian. Thus we re-scale the momentum and the 
Hamilton operators as p± — 2v p± /71 , and H = 2H/ji , and then the Hamiltonian takes the following form (we omit 
the tilde): 



H 



/o 








P-\ 








p+ 








p- 





2 


w 





2 


0/ 



(42) 



To calculate the time dependence of the operator x(i) we shall use Eq. After some matrix algebra we found 
the following result for the decomposition of the Hamilton operator (|42p into a sum of projection operators: H = 
Y.t^iEaQa, where 



El 



Qi 



I A — R Ia 



T 



Q2 



h + R h-T 



L-R h-T 



^Ji + P+P- = VT+p 

h + R h + 



Q4 



T 



(43) 
(44) 



R = 



I e-^'-^ \ 

e^^'P 

1 

V 10/ 



T = 










p- 




1 





-1 









p+ 





1 







\ 


p- 





1 / 



and p± ~ IpI e^"''. 



(45) 



and /4 is the 4x4 unit matrix. One can easily show that R^ — = I4, ie, both R and T are involutory operators. 
Similarly, R and T commute with each other, ie, [R, T] = 0, and the four operators Qi, . . . , Qi are indeed projection 
operators satisfying the usual relations: QaQb = Sat Qa and X]a=i ~ -^4- One can also show by direct calculation 
that H = Ea Qa holds. 

We now can use Eq. Q or to find time dependence of the position operator. Note that to obtain the explicit 
form of the Zittcrbewegung amplitudes it is more useful to apply Eq. ([5|) since it is trivial to calculate the derivative 
of the Hamiltonian (|42|) with respect to the momentum operator p. After a tedious but straightforward calculation 
we have 



x(t) = x(0) + Wt + Z(i), where 



EdEi dfl p p 



dp 



dp 





J 

P+ 



P 



p_ 


1 





(46) 



(47) 
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and 



h e X p 



{cos{2 at) - 1) 



13' 





-pp- 













13 P+ 


PP+ 













pp- 





-P^ 



+ (cos(2/3t) - 1) 



+ (cos(2f)-l) 



a2 













-a^ 





— ap 







p2 








— ap_ 


. 


-P^ 


P^ 





2p- 








-2p^ 





-2p+ 



isva{2j3t) 







/ 







i sin(2 


at) 
























-0? 


pi 





—ap- 









ap+ 








— Q:p_ 





-P^ 


ap+ 







P^ 






-^P- 



pp+ 



^2 



2p+ 






-2p. 



-2p2 





2p2 



+ i sin(2 1) 



f 


^2^+ 



-2p2 


-2p- 






2p+ 


-2p2 



-2p_ \ 


2p^ 




V 









Q.p- 
















Q.P+ 














-Q.p- 









+ i sin(2f2t) 



/ 






-P- 




/ 



-p+ 



^2 



-p- 



p' 







p' 





(48) 



where e = (0, 0, 1), a = — 1, /3 = f2 + 1 and t is in units of h/^\. One can easily see that the above operator x(t) 
is a hermitian operator. The vector e x p = (— Py,Pa;,0) is perpendicular to the momentum vector p = {px,Py,Q)- 
Therefore, the operator x(t) has one longitudinal and three transversal modes parallel and pcrpcndiclar to the mo- 
mentum p, respectively. It is clear from the result that the oscillatory term Z(<:) corresponding to Zitterbewegung is 
a sum of individual oscillatory terms with four different frequencies. 



